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We investigate the addition spectrum of disordered quantum dots contain- 
ing spinless interacting fermions using the self-consistent Hartree-Fock ap- 
proximation. We concentrate on the regime r^ ^ 1, with finite dimension- 
less conductance g. We find that in this approximation the peak spacing 
fluctuations do not scale with the mean single particle level spacing for ei- 
ther Coulomb or nearest neighbour interactions when r^ ^ 1. We also 
show that Koopmans' approximation to the addition spectrum can lead 
to errors that are of order the mean level spacing or larger, both in the 
mean addition spectrum peak spacings, and in the peak spacing fluctuations. 



I. INTRODUCTION 

We consider the response of the ground state of spinless fermions to the addition or 
removal of a particle. To this end, we apply the self-consistent Hartree-Fock (SCHF) ap- 
proximation: a non-perturbative effective single particle theory. 
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Koopmans' theorem |]T| states that the single particle SCHF energy levels describe the 
affinity and ionisation energy spectra for the unoccupied and occupied states respectively. 
The approximation involved is that all the other particles do not react to this process. This 
approximation is generally considered to be good when the single particle wavefunctions are 
extended: corrections to each wavefunction due the rearrangement of the system following 
the addition or removal of a particle are expected to be 0{1/N), where A^ is the number of 
particles in the system 0. Moreover, these corrections to the wavefunctions are expected 
to disappear in the limit of vanishing disorder (having e.g. periodic boundary conditions), 
where the single particle wavefunctions are free waves, as well as in the limit of vanishing 
interaction. It is then loosely assumed that in the thermodynamic limit, Koopmans' theorem 
becomes exact even for disordered systems, and should be sufficiently accurate for mesoscopic 
samples. It is evident however, that if the physical quantities at hand require an energy 
resolution of order the mean single particle level spacing. A, the validity of Koopmans' 
theorem should be reconsidered. 

Our analysis of Koopmans' theorem for a quantum dot is closely related to the addition 
spectrum of the latter: the spectrum of energy differences between states with total particle 
number different by unity. We consider only the energy differences between ground states, 
which is experimentally accessible through resonant tunnelling P-^, and capacitance P,|^ 
measurements at low bias and temperature. We stress that while our analysis here pertains 
to some aspects of the experiments, a direct comparison is not feasible: ffistly we consider 
spinless electrons, and secondly we consider disordered systems in the diffusive regime; in 
Ref. the mean free path is of order the sample size, whereas in Ref. [Q the mean free 
path is much larger than the system size, and ergodicity is ensured by a chaotic boundary 
shape. 

The position of the observed resonant tunnelling (RT) conductance peaks can be related 
to the ground state energy difference ^n = Eg{N, Vg) — Eg{N — 1, Vg) where Eg{N, Vg) 
is the ground state energy of the dot with N electrons, at gate voltage Vg- The spacing 
between consecutive peaks is thus related to 



A2(iV) = Eg{N + 1, Vg) - Eg{N, Vg) - Eg{N, V^) + Eg{N -1,V^) . (1) 

Within the constant interaction (CI) model the ground state energy is simply the sum of 
filled single particle energies e(n) plus A^(A^ — l)Vo/2, where Vq is the constant interaction. 
Taking Vg = Vq, the peak spacing trivially reduces to 

A2{N)=e{N + l)-e{N) + Vo (2) 

and so, in the diffusive regime, displays shifted Wigner-Dyson (WD) statistics HJ^ up to 
corrections in one over the dimensionless conductance, g ^^: P{s) = 7rs/2exp(— 7rs^/4) for 
zero magnetic field, the case that we consider here; s = (A2 — Vo)/A. 

Recent experiments on quantum dots f^-^l have shown that whilst the mean peak spac- 
ings are well described by the CI model, the fiuctuations are not described by Wigner-Dyson 
statistics. It is found that the distribution of A2 is roughly Gaussian 0,0, with broader 
non-Gaussian tails seen in Ref. ^. In Refs. [^Q the variance of the fluctuations was found 
to be considerably larger than that given by the WD distribution. Further experimen- 
tal observations, including correlations of peak heights ITTJ-^ and the sensitivity to an 
Aharanov-Bohm flux |T^JTI|JT^ are not consistent with random matrix results, and suggest 
a breakdown of the naive single particle picture. To investigate this, one needs information 
on the ground state wavefunction. 



Blanter et al [T^ have evaluated the fluctuations within a Hartree-Fock framework, ne- 
glecting effects due to the change in gate voltage Vg to Vq. To this end they applied the 
Random Phase Approximation (RPA) to generate the screened interaction in the conflned 
geometry, and assumed that all HF level spacings, except for the Coulomb gap, are de- 
scribed by WD statistics. Implicitly assuming Koopmans' theorem [0 to be valid, and using 
wavefunction statistics established for non-interacting electrons in a random potential, they 
calculated the fluctuations of A2 beyond the CI model. These additional fluctuations were 
found to be parametrically small (in 1/g) and proportional to A. Hence, the total fluctu- 
ations in A2 were found to be proportional to A. The analysis of Ref. []IB[ is consistent 



in the limits g ^ 1, Vg <^ 1. The parameter r^ characterises the relative importance of 



interactions in the electronic system, and is defined as the mean electron separation in units 
of the effective Bohr radius. 

Exact numerical calculations on small disordered dots [Q did produce large Gaussian 
distributed fluctuations at experimental densities. It was claimed that for strong enough 
interaction 5A2/ {/S.2) is universal, independent of the interaction strength and disorder, 
where 5/S.2 denotes the typical (RMS) size of the fluctuations, and the angle brackets denote 
disorder averaging. This universal constant was found to be approximately .10 — .17, to be 
compared with the WD result for the CI model: .52A/(A + Vq). We note that the typical 
experimental value for the charging energy Vq, is much larger than A. The scaling with (A2) 
suggested by this analysis is in stark contrast to the scaling with A obtained in Ref. [jl6[. 

Stopa has considered ballistic chaotic billiards numerically, using local density functional 
theory |T^. In this case, it was claimed that the fluctuations arise due to strongly scarred 



wavefunctions in the self consistent potential. As a result of these scars, an asymmetric 
distribution of A2(A^) was found, including strong correlations over A^. It was then further 
noted that what is actually measured (i.e. the change in the gate voltage between resonant 
tunnelling peaks) is not simply related to A2(iV) when the dependence of A2(A^) on the gate 
voltage is strong. It was then claimed that a self-consistent calculation of AVg with A2(A^) 
retrieves a symmetric distribution of A2(A^), and reduces peak to peak correlations. 

A further suggestion that the coupling to the gate is important in understanding the 
fluctuations has been made with reference to the CI model, with WD statistics for the single 
particle levels |T^ . The authors claim that the required distribution of A2 can be generated. 



except for the non-Gaussian tails, through the de-correlation of neighbouring levels under 
a parametric change in the Hamiltonian (mediated by Vq)- However, the degree of de- 
correlation induced by AV^ is left as a fitting parameter. 

In this paper we present numerical calculations within the SCHF approximation, con- 
sidering larger samples than is feasible by exact diagonalisation 0. This approximation 
has been seen to be quite good for the calculation of persistent currents in similar systems 
il9| . We show that fluctuations large compared to the single particle level spacing can arise 



without recourse to varying the sample shape, size or gate to dot couphng, supposing these 
to be additional effects. We further demonstrate that approximating the addition spectrum 
spacings by applying Koopmans' theorem can lead to large errors in the calculation of the 
spacing statistics. 

We consider separately both a long range (Coulomb) bare interaction and a short range 
(nearest neighbour) bare interaction. In section || we introduce our model in detail; in 
section |T| we present a short discussion of the implications of Koopmans' theorem; in 



section |^ we present and discuss our numerical results, which are then summarised in the 



final section. 

II. THE MODEL 

We address the following tight binding Hamiltonian for spinless fermions 

H = J2 Wicfci - t Y. ct+^c^ + Y E Mijijctc^Cjd (3) 

i i,ri ij 

where i is the site index, t] describes the set of nearest neighbours, Wi is the random on site 
energy in the range [—W/2,W/2], and t the hopping matrix element, henceforth taken as 
unity. We study separately, both a Coulomb interaction potential, 

M,,,, = l/|r,-r,| (4) 

and a short range potential plus a constant term M^ (see below), 

Mijij = (5,,i+^ + MJ . (5) 

We consider a 2D system with periodic boundary conditions, and choose to define 

\ri - Tjl^ = [Llsm^{7in^/L^) + Llsm^{7iny/Ly)]/n'^ , (6) 

where {rix, Uy) = ri — Tj. 

All distances are measured in units of the lattice spacing a; the physical parameters are 
therefore Uq = e^/a, t = h^ /2ma'^. The standard definition for r^ is given, for low filling, 



by Tg = Uq / {t^/Am}) , where v = N/A is the fiUing factor on the tight binding lattice with 
A sites. The dimensionless conductance g can be approximated, again for low filling, using 
the Born approximation. We find g = 967Tu{t/W)'^, which is valid for A,N^g^ 1. Here 
u ^ 1/4 throughout, so that r^ ^ .56f/o/t, and g ^ 75{t/W)'^. 

Having identified the parameters of our model with the standard ones employed in the 
theory of a continuous electron gas, we note that in the limit of small r^ and 1/g the leading 



order term for the typical interaction dependent fluctuations predicted by Blanter et al ||T^ 
is ~ UoA/ty/g. For the torus geometry considered here, this contribution, being a surface 
term, vanishes identically. Their prediction then reduces to typical fluctuations in addition 
to those of the CI model to be of order UoA/tg. 

The torus geometry has the advantage over geometries with hard walls whereby in the 
former, the compensating background charge provides a trivial shift in all the site energies, 
and can be removed. In a bounded dot, with an overall charge, the excess charge may build 
up near the boundary, depending on the position of nearby metallic plates and gates. These 
effects are geometry specific |T^ . Upon adding an electron, the average charge configuration 
may change considerably (the configuration is strongly geometry dependent). As the gate 
voltage is varied to allow for the next electron addition, the background potential could 
have changed causing further charge rearrangement. Whilst it is of great interest to analyse 
this issue (which may play an important role in the peak spacing fluctuations as well as 
undermining the naive single particle picture by further reducing the accuracy of Koopmans' 
theorem ljl3|); we concentrate here on effects due entirely to the intrinsic rearrangement of 
the dot. From this point of view our analysis may be taken as an attempt to establish an 
upper bound criterion for the breakdown of Koopmans' theorem. In reality it may break 
down earlier due to other non-universal factors. During the completion of this work very 
recent experimental evidence for significant rearrangement has been produced |T^. It is 
argued that rearrangements due to adding an electron are far greater than rearrangements 
due merely to a change in shape. 



When considering the short range interaction, the mean charging energy Vq in Eq.(^ 
must be put in by hand through M^. of equation (^). The way in which this is done depends 
on the physical situation being modelled, and is highly geometry dependent (vis-a-vis the 
gates). We stress that the value of Mc does not affect the physical results. We choose to 
insert 

— , r — r' 

r ,r' I I 

This value for Mc is defined such that if the charge is uniformly spread over the dot, the 
average charging energies in the Coulomb and nearest-neighbour cases roughly coincide. This 
choice has been made for simplicity, but corresponds to the premise that the interactions of 
the iV-electron gas with the positive background and with itself is the same for both models 
considered. Exchange contributions, which tend to reduce the total charging energy, are 
included insofar as to cancel both the on-site contributions to the energy, and the unphysical 



self- interaction of electrons, but are otherwise neglected [T^JS^ . The energy associated with 
charging the system uniformly is UqN{N — 1)/{2A'^) J2'^ ^, |r — r'|^^ in the Coulomb case, and 
UoN{N — 1)/{2A'^){AA + McA"^) in the nearest-neighbour case. Equation (0) follows from 
equating these energies. This estimate can be systematically improved if the above premise 
is taken as the definition, not only by correctly accounting for the exchange contributions, 
but also by considering single particle wavefunction statistics in the diffusive regime. In this 
case wavefunction correlation functions such as (|?/'j(r)p|^j(r')p) are required in order the 
calculate the average electrostatic energy, where here and after (...) denotes averaging over 
the disorder ensemble. 



In Ref. [0 it was assumed that the (RPA) screening can be taken into account before 
constructing the Slater determinant ground state, and therefore their result corresponds to 
a short-ranged effective interaction. It is not clear that this remains a consistent procedure 
when calculating the ground state energy self-consistently. The reason for the inconsistency 
is that many of the diagrams generated by the SCHF approximation are already included 
in the RPA calculation of the screening, resulting in double counting. On the other hand, if 



the screening is generated externally (e.g. by close metallic gates), then it is consistent to 
insert a short ranged bare interaction, and this is the point of view taken here. 

In some sense, the Coulomb interaction results can be considered as the opposite limit of 
the nearest-neighbour interaction, and is of interest in this context. However, it is more diffi- 
cult to physically motivate the use of a Coulombic bare interaction unless one is considering 
very low electron densities. Screening is indeed weak in a 2d electron gas in a vacuum, even 
at high density, but the SCHF procedure cannot correctly generate screening by itself; whilst 
it can screen the Hartree contributions (as discussed above), it does not screen the exchange 
(Fock) term. However, we have verified that for the range of parameters considered here, 
fiuctuations in the Hartree energy are larger than the typical fiuctuations of the exchange 
energy. This suggests that the error made in not screening the exchange term correctly is 
not overly important. 

III. IMPLICATIONS OF KOOPMANS' THEOREM 

Let us now consider the form of A2, and approximations to it given by applying Koop- 
mans' theorem. We denote the diagonal matrix elements of the one-body operators in (|^) by 
Tj^, and the antisymmetrised Hartree- Fock interaction by V^ , where hereafter the subscripts 
denote single particle states, and the superscript A^ denotes the number of particles present 
and identifies the self-consistent basis of single particle wavefunctions being employed, ipf . 
For the torus geometry, where the gate voltage and background potential represent a trivial 
shift that can be omitted, the SCHF ground state energy is given by 

where ej" is the /th SCHF single particle energy for a system of m particles in the ground 

state: 

m 



Using §), we find (c.f. Eq.(|)) 

N 

+ E (^#i^- - ^#7') + ^ E (^r ' - 2^.f + yt') ■ (10) 

Applying Koopmans' approximation corresponds to dropping the superscripts and employ- 
ing an appropriate fixed basis. The theorem implies that the effective single particle states 
do not depend on the occupation of these states. In particular, Koopmans' theorem yields 
ejv+i fo^ ^^6 minimum energy required to add a particle to a system of A^ particles, and e^ 
for the maximum energy gained by removing a particle from the same system; in both cases 
the final state is a ground state. Clearly ef^ as well as the ground state energy depend on 
m, even in Koopmans' approximation, through the number of terms in the sum in Eqs.(^ 
and (|^) respectively. It is then easy to see that Koopmans' approximation yields 

A',^{N)=e^_,,-e^. (11) 

We also consider two other approximations to A2 that involve calculating two self-consistent 
bases rather than just one: 

Ar(iV) = 4+i-4-' (12) 

A^^iV) = 4:1 - 4 • (13) 



All three estimates (p!l|-|T3|) coincide with A2(A^) of Eq.(|TOD if Koopmans' theorem holds. To 



connect with the notation of Ref . [0 , and to demonstrate the difference between the above 



three approximations and the fully self-consistent result, we provide a schematic diagram of 
the SCHF spectra in Fig. |I|. 



'AAN+D 



A,(N) AJN) 



N N+1 

FIG. 1. A schematic digram of the SCHF spectra of N and A^ + 1 particles. 



Since the self-consistent basis of A^ particles provides the lowest energy for A^ occupied 
levels, and similarly for A^ — 1 particles, the following relations are clear: 

N-l 



N-1 1 N-l N-l 1 

i ij i 



N~l 
2 ^ • *j 



U 



N 1 TV N 1 Af 

i ij i ij 



(14) 



Combining these equations, we find that Ae(A^) = A2^(A^) — A2^(A^) > 0, or equivalently 



Ae(A^) 



^7V-l _ N 

-N ^N 



eZ > . 



(15) 



The equalities in (|14D,(|T5D only hold when no modification of the effective single particle 
wavefunctions occurs following the addition of an electron. In a disordered dot, in which 
there are no spatial symmetries, such a modification will always take place, and so Ae can 
be considered strictly positive. 



N+1 
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FIG. 2. A schematic diagram representing the expectation value of the Hamiltonian in the spaces 
of A^ — 1, N and A^ + 1 Slater determinants (superimposed). Distances are only meaningful within 
a given space. In order to define distances between two Slater determinants in two different spaces, 



i.e. ^^ and vj/^+i^ we include the first unoccupied state with ^^ |21|. In this way Koopmans' 
theorem is exact when the two Slater determinants coincide. The SCHF solutions correspond to 
minima in these surfaces. It can be seen that the Koopmans' approximations e]^_^_i and e}^ to the 
addition energy Eg(N + 1) — Eg{N) are upper and lower bounds respectively. Similarly for the 
addition energy Eg{N) — Eg{N — 1). The definition of the Koopmans' approximation to A2 given 
by Eq.(^) can be seen to contain the difference between these bounds. 

The difference Ae provides a measure of the effectiveness of Koopmans' theorem. To 
demonstrate this we present, in Fig.^, a schematic diagram of the surface of expectation 
values of the many-body Hamiltonian in the space of Slater determinants of A^ — 1, A^, and 
A^ + 1 particles. The SCHF ground states correspond to minima in these surfaces. From 
the diagram, it is clear that the energies ej^_,_^ and e^^^ are upper bounds to the respective 
addition energies, and the energies t%~^ and e^ are lower bounds to the addition energies. 
The approximation A2S is therefore obtained by subtracting a lower bound (ej^) from an 
upper bound (ctv+i)- ^s a result, the average value contains the average difference between 
the two bounds in addition to the correct mean A2. It is generally assumed that the difference 
between the two bounds vanishes in the thermodynamic limit, and therefore so does Ae. We 
shall see that our results do not show any indication that this is the case. On the other hand, 
A2^ corresponds to the difference of two upper bounds to the two relevant addition energies. 
Regardless of the quality of the upper bound, so long as it is not strongly dependent on 
the number of particles present, both the particle number and disorder averaged results are 
good. The third approximation to A2^ corresponds to the difference of two lower bounds, 
and like A2^ is good in the mean. It is for this reason that we introduce these alternative 
approximations. It is easy to see that A2^(A^) — A2^(A^) = Ae(A^ + l) and therefore provides 
no further information. On the other hand, the fluctuations of A2^ can be different from 
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fc2 



those of A2^ , and so are investigated separately. We note that in a clean system at r^ below 



the Wigner crystal transition [^], the minima would align in Fig. 0, reflecting the validity 
of Koopmans' theorem in that limit. 

Let us briefly discuss the non-self-consistent single particle picture, for which the Koop- 
mans' approximations ([TTl - |l3l) and Eq. ([lO|) all coincide: 



7V-1 

A2(iV) = T^+i - T^ + E {Vn+ij - V^,) + V^+i^ . (16) 

j 

Here, the term non-self-consistent approximation refers to a scheme where a set of effective 

single-particle states is given (e.g. by solving the N-electron SCHF problem), and utilised 

for any number of particles present in the system. The nearest neighbour spacings between 

levels that are both occupied or unoccupied has a similar form: 

N 
^m+1 ~ ^m = Tm+1 — T^ + 2_^ (Ki+lj " Knjj , (17) 



the major difference between (|T^) and ( |T7| ) is the additional unbalanced matrix element 
Vat+iat appearing in ([T6|) . Let us also suppose that in this simple single particle scheme 
the electrons interact with a short-ranged effective interaction. Blanter et al ^^ introduce 



the hypothesis that the (normalised) spacings (|T^ and A2 — Vat+i jv obey WD statistics up 
to corrections in 1/g. Further assuming that the wavefunction correlations are still close 
to those of non-interacting particles leads, for the short-ranged effective interaction, to the 
result Var(Vij) ~ (UoA/tg)^ [^], so that the interaction dependent contribution to SA2 
scales like UoA/tg [jl6[. This analysis is valid in the regime r^ <^ 1 and g ^ 1, implying 
that Koopmans' theorem is a good approximation in that regime. 

IV. RESULTS AND DISCUSSION 

In this section we present and discuss the results of the numerical simulations for both the 
nearest-neighbour and the Coulomb bare potentials. To make each subsection self-contained 
there is some repetition. 
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A. Short Range Interactions 



We consider first the case of a nearest neiglibour bare interaction potential as defined in 
Eq.(^. We begin by plotting the distributions of both the level spacings ([T7[) and the gap 
(p^) of the SCHF spectrum for finite r^, g. In this case (16) and (|l^) are calculated in the 
self-consistent basis of A^ particles. In Fig. |^ it is seen that the normalised level spacings 
between occupied states show an increasing deviation from WD to Poisson statistics as Uq 
is increased. This is also true for the unoccupied states, but to a much greater extent. 
The difference between occupied and unoccupied states in the SCHF approximation will be 
discussed in greater detail later. We interpret the tendency towards Poisson statistics as a 
signature of the incipient localisation of the effective one particle states. 

The normalised gap ( Ag^ ) distribution tends towards a more symmetric distribution that 
is approximately Gaussian as Uq is increased. 
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FIG. 3. SCHF level spacing distributions P[s] for a) occupied states, b) unoccupied states; 
s = AE/{AE), and c) the gap AgS where s = {A^^ - {A2^))/5A\^. The solid lines show the WD 
distribution, and in (c) the dashed line follows a Gaussian law. The samples were 8*9 lattices 
with 14 electrons and nearest-neighbour interactions; W = 2. r^ ~ 0.56C/o/t. The statistics were 
obtained from an ensemble of 2500 samples. 

We have also investigated the gap (A2) distribution obtained within the fully self- 
consistent scheme, which we show in Fig. ^. We find that as Uq is increased, the distribution 
evolves from a WD form to a more symmetric distribution similar to a Gaussian. 




FIG. 4. Distributions P[s\ where s = (A2 — {A2))/5A2 for various interaction strengths obtained 
self consistently for the nearest neighbour interaction. The solid line show the WD distribution, 
the dashed line shows the Gaussian distribution. The samples were 7*8 lattices with N = 15, 
W = 4, and the statistics were obtained from an ensemble of 3000 samples, r^ ~ 0.56Uo/t. 

We shall concentrate first on the mean values of these distributions. A typical dependence 
of (A2) on the interaction parameter Uq is plotted in figure]^. Whilst (A2^) and (A2^) provide 
a good approximation, we see a strong deviation of (A2^). For all the system sizes considered, 
this effect occurs at r^ ~ C'(l). Results for the CI model, evaluated as described above, 
are plotted for comparison. Deviations of order 0{Uo/A) from the CI model appear above 
f/o ^ 2 (r, ^ 1). 

14 



Elsewhere [^ , we show that the ground state develops large density modulations as Uq is 



increased beyond r^ ~ C>{1). These ground state charge density modulations (CDMs) explain 
the deviations of (A2) from the CI model prediction: they reduce the average addition energy 
by up to 4Uo/At when u < 1/2 for a commensurate lattice. When u > 1/2 the mean charging 
energy can be correspondingly increased. 

We are also now in a position to understand why the level spacing statistics between 
unoccupied states show an increased tendency towards a Poisson distribution: the density 
modulations that appear in the ground state alter the potential felt by the unoccupied 
states. These modulations are not spatially ordered, as is demonstrated in Ref. p4|. Hence, 



as Uq is increased, the unoccupied states see an effective potential with increasingly strong 
modulations and tend to localise, whence the tendency towards a Poissonian distribution. 
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FIG. 5. Typical result for the mean Coulomb gap averaged over 400 disorder realisations. Here, 

W = 4, the lattice is 9 * 8 and iV = 15. The dashed line is the CI result, r^ « 0.62Uo/t. 

Let us consider the error in (A2^) in more detail. As can be seen in Fig. ^, (Ae)/A 
increases with the system size for lattices up to about 7*8; for larger systems {A ^ 50) it 
seems that the error becomes proportional to A, and is of order A when r^ is of order unity. 
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FIG. 6. (Ae)/A against the sample area A after 75 to 1000 disorder realisations. The dotted 
lines are proportional to A as guides for the eye. r^ ~ Uo/2t. 

We find that the nature of the disorder dependence of (Ae) depends on the interaction 
strength, as seen in Fig. |^. The change in dependence occurs at interactions strengths 
corresponding to r^ ^ 1 for all the sample sizes considered. One might be surprised that 
deviations from Koopmans' theorem do not smoothly decrease as W —>■ since, in the 
limit of vanishing disorder, Koopmans' theorem becomes exact on the torus due to the 
restoration of translational symmetry. However, when VT — > 0, the spectrum develops many 
near degeneracies such that the effective perturbation due to Uq is magnified as W -^ 0. 
In the limit r^ — ^ 0, f^ ^ 1, the typical size of the matrix elements Vijki which drive the 
rearrangement scale like SVijki ~ ^sA/g [^, thus one expects that in this regime (Ae) 
should increase with disorder. We find only a weak increase with disorder for r^ ^ 1. 

In figure ^ we plot the interaction dependence of (Ae)/A. We find that at small Uq, 
(Ae)/A oc (f/o/t)^, with deviations for larger Uq. In fact, this quadratic behaviour can 
be understood using second order perturbation theory. To see this we refer back to the 
schematic diagram of Fig. §. A shift |^ occurs in the ground state configuration when 
a particle is added, which is represented by a misalignment of the minima. This shift is, 
to leading order, linear in Uq. Since the SCHF ground state energy is a minimum in the 
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expectation value of the Hamiltonian, the difference of ground state energies for the two 
configurations will be quadratic in this shift. Furthermore, the local curvature tensor is 
independent of Uq when the interaction matrix elements are small compared to the mean 
level spacing [^. Thus, (Ae) scales like Uq in the perturbative regime. The indication is 
that second order perturbation theory is qualitatively good even for r^ ~ 1. We note that 
since both the shift and the local curvature tensor depend on the disorder, there is no such 
simple W dependence. 
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FIG. 7. (Ae)/t against disorder W averaged over 200 disorder realisations for a 11 * 10 lattice 
with N = 28. rs « OMUo/t. 
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FIG. 8. (Ae)/A against Uo after 75 to 1000 disorder realisations, for a range of sample sizes: 
u ss 1/4, on a log-log scale. The dashed line is a plot of (Ae)/A oc 11^/1^, r^ ss Uo/2t. 

To summarise the results for the mean Coulomb gap, we find that Koopmans' approxima- 
tion (|Il|) makes an error in the mean charging energy which for small Vg, A '^ 50, and fixed 
disorder scales like (Ae) oc r^. There is also evidence that for the larger sizes [A > 50), far 
beyond that accessible by exact diagonalisation, that (Ae) oc r^A. The latter dependence 
is consistent with the expectation that for sufficiently small r^, 1/g, perturbation theory 
is valid when the effective interaction is short-ranged. We find that (Ae) ~ 0{A) when 



r^ ~ (9(1). To understand this result, we return to (|T^), i\l7\} , and fix the basis to be the 
self-consistent one for A^ particles, so that (p!6|) now describes A^^. Since we have verified 
that the level spacings (|1^) show nearest neighbour separation statistics which are close to 
WD for all m ^ N, with an approximately constant density of states, we are led to conclude 
that (Ae) arises due to the fundamental difference between occupied and unoccupied levels 
in the SCHF approximation. In short, whilst {V^_^_ij — V^j) for m ^ N vanishes as expected, 
(V^Ai, - ^Nj) does not. Indeed (Ae) oc Er'(^Ai, - V^,)- 
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FIG. 9. Typical result for the typical fluctuations of the Coulomb gap averaged over 400 disorder 

realisations. Here, VF = 4, the lattice is 9 * 8 and N = 15. r^ ~ 0.62f^[)/t. 

Turning now to the fiuctuations in A2, we plot an example result in Fig. p. The most 



striking behaviour is the asymptotic saturation of the fluctuations (verified but not shown 
for even stronger interactions). As with the deviations of (A2) from the CI model, this oc- 
curs over the same range of interactions for all the system sizes considered, and is associated 
with the appearance of CDMs. Over the range of interaction strengths shown the fluctua- 
tions have not completely saturated, but the ground state density modulations are already 



present |]2J], such that, at low filling, the short-ranged contribution to the interaction energy 
is reduced. In the limit of strong interactions the charge segregates at a kinetic energy cost 
of order 0(t), and Uq plays no further role in the ground state energy fluctuations. The 
fluctuations therefore become sub-linear in Uo, and eventually saturate to an interaction 
independent value. Moreover, the results for the fluctuations become strongly geometry and 



filling factor dependent p6|. We note that the observed saturation is in fact an artifact of 



the sharp cut-off in the interaction range: with a longer-ranged interaction, charge segre- 
gation cannot eliminate contributions due to the interaction (although it may significantly 
reduce them), and the fluctuations would no longer be bounded simply by kinetic energy 
considerations. 

In Fig. y it can also be seen that for strong interactions, the fluctuations are overesti- 
mated by A2^ and Ag^, and underestimated by A2^. This can be understood within the 
picture given above of charge density modulations. In this case, an occupied state that is 
removed non-self-consistently will yield less energy than can be gained when the system is 
allowed to reorganise, but the typical size of this error saturates to an interaction indepen- 
dent value for the same reason that the SCHF fluctuations do. If on the other hand an 
unoccupied state is occupied non-self-consistently, it is not possible to avoid contributions 
from the short-ranged part of the potential (we do not consider strongly Anderson localised 
states at very low filling), and the typical error increases indefinitely with the interaction 
strength. As a result Ag^ underestimates the fluctuations by an amount that saturates to an 
interaction independent value, whereas fluctuations in the charging energy predicted by A2^ 
and Aj^ grow with Uq indefinitely; the errors made in employing the latter approximations 
diverge with the interaction strength. 
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FIG. 10. 5IS.2/ IS. against A, with W = A and v ~ 1/4, after 75 to 1000 disorder realisations. The 
RMT results is shown as a dotted line, the dashed line is proportional to L and the dot-dashed 
line is proportional to A. rg ~ Uo/2t. 

Concentrating now on the fully self consistent results, the fluctuations in the charg- 



ing energy are plotted against the sample size in figure |10|. As expected, for very weak 
interactions, the typical fluctuations vanish like 1/A, being dominated by kinetic energy 
fluctuations. For stronger interactions, this dependence no longer holds: for r^ ^ 1 our 
results are in broad agreement with Ref. ||16[, but do not agree with their suggestion that 
the typical fluctuations remain proportional to A for r^ > 0{1). This appears to conflict 
with a simple single-parameter scaling argument |^ . The appearance of fluctuations that 
do not scale with A coincides with the appearance of density modulations. We stress that 
in this model there can be no physical connection between the amplitude of the constant 
interaction and the amplitude of the fluctuations. 

For strong interactions the dominant disorder dependence appears to develop only for 
W ^ 4, where it is consistent with the emergence of a linear dependence to be expected 
from spatial rearrangements in the disorder potential. An example is plotted in Fig. |ll]. We 
reiterate that we also find strong geometry and filling factor dependences. It is extremely 
difficult to extract disorder scalings in such small systems because W/t is required to be fairly 
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large to generate diffusive motion, which in turn stretches the spectrum in the tails. This 
can be seen at weak interaction, where one would have hoped to see a disorder independent 
plateau in A (i.e 6A2 at Uq = 0). 
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FIG. 11. A typical result for 5/S.2/t against W ^ averaged over 200 disorder realisations. Here the 
sample is 11 * 10, iV = 28.rs « 0.54C/o/t. 

To summarise then, we find 5A2 ~ 0.52 A + ar^A + O(r^), where a is an undetermined 
constant or function of disorder strength. We note that the disorder scaling is not clear 
because of the residual dependence of A on W^. 



B. Long Range Interactions 



We consider here the results for the Coulombic bare potential. 
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FIG. 12. SCHF level spacing distributions P[s] for a) occupied states, b) unoccupied states; 
s = AE/{AE), and c) the gap AgS where s = {A^^ - {A2^))/5A2^. The solid lines show the WD 
distribution, and in (c) the dashed line follows a Gaussian law. The samples were 8*9 lattices 
with 14 electrons and Coulomb interactions; W = 2. r^ ~ 0.56f/o/t. The statistics were obtained 
from an ensemble of 2500 samples. 

We first study the distributions of both the level spacings ([l^) and the gap ( |16|) of the 
SCHF spectrum at finite r^, g. In this case ([T6| ) and (|l3) are calculated in the self-consistent 



basis of A^ particles. In Fig. |T2| it is seen that the normalised level spacings between 
occupied states obey statistics very close to WD for all interaction strengths considered. 



Between occupied states (Fig. ^a) there is a mild deviation towards Poisson statistics 
for the strongest interaction strengths, indicative of a weak tendency towards localisation. 
Between unoccupied states (Fig. |12|b) the distribution is even closer to WD. On the other 
hand, the normalised gap distribution clearly tends towards a more symmetric distribution 
that is approximately Gaussian. 

We have also investigated the gap (A2) distribution obtained within the fully self- 
consistent scheme, which we show in Fig. |T3l Again we find that as Uq is increased, the 
distribution evolves from a WD form to a symmetric distribution similar to a Gaussian. 
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FIG. 13. Distributions P[s] where s = (A2 — (A2))/(5A2 for various interaction strengths obtained 
self consistently for the Coulomb interaction. The solid line show the WD distribution, the dashed 
line shows the Gaussian distribution. The samples were 7*8 lattices with N = 15, W = A, and 
the statistics were obtained from an ensemble of 3000 samples, r^ ~ 0.56C/o/^- 



We shall concentrate first on the mean values obtained from these distributions, and 



will come to the variance later in the section. Figure |14| shows a comparison of (A2), with 
the various approximations to it, plotted against Uq. Whilst {^2^) ^iid (A2^) provide a 
good approximation to (A2), we see a clear deviation of (A2^). Results for the CI model, 
evaluated as described above, are plotted for comparison. That the CI model is good in the 
mean indicates that the single particle wavefunctions remain roughly uniformly distributed 
over the dot for all r^ considered. 

In figure |T5| we plot (Ae)/(A2) against the sample area A, for an intermediate disorder 
strength {W = 4). Since we know from Fig.[T^ that (Ag^) ~ (A2), then (Ae) is very close to 
(Ag^) — (A2), the total error made by applying Koopmans' theorem. For A ^ 50, we find that 
for a fixed interaction strength (Ae) oc L ((Ae)/(A2) oc A). For larger samples with Uq ^ 2 
we find a weakening in the dependence, but see no indication that it will vanish relative 
to A. The result that the deviations from Koopmans' approximation increase with system 
size (when compared to A), showing no sign of saturation, is admittedly strange, and may 
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be an artifact of the specific model considered here. However, the result that Koopmans' 
approximation appears to fail even as the system size tends towards the thermodynamic 



limit, is in line with our findings for the short-ranged case p8[] 
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FIG. 14. Typical results for the mean Coulomb gap under various approximation schemes, aver- 
aged over 400 disorder realisations. Here, W = 4, the lattice is 11 * 10 and N = 28. r^ ~ O.bAUo/t. 
The dashed line is the CI result for the mean. 




FIG. 15. (Ae)/(A2) against the sample area A] W = A, v k, 1/4, and the results averaged 
over 300 to 1000 disorder realisations, for a range of interaction strengths. The dashed lines show 
(Ae)/(A2) oc j4 as guides for the eye. r^ « JJ^jlt. 
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It is interesting to see how the error depends on disorder. In figure we plot (Ae)/t 
for a range of disorder strengths: the disorder dependence as a function of interaction is 
weak, but not simple. Similarly to the short-ranged case, the deviations from Koopmans' 
approximation do not decrease for small disorder. This occurs for the same reasons as for 
the short-ranged case. Here too the typical size of the matrix elements Vijki which drive 
the rearrangement scale inversely with g. One thus expects that in this regime (Ae) should 
increase with disorder, and this is indeed seen in the figure. For Uq ^ 4t, Ae decreases with 
disorder at sufficiently large W, with evidence of a turning point {dAe/dW = 0) at Uq = 4:t, 
W ^At. 
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FIG. 16. (Ae)/A against disorder W/t averaged over a ensembles of 300, 11 * 10 samples with 
28 particles, r^ ~ O.BAUo/t. The dashed line is proportional to yW. 



In figure |T^ we plot the interaction dependence of (Ae)/A. We find that at small Uq, 
(Ae)/A oc (f/o/t)^, with deviations for larger Uq- This quadratic behaviour has the same 
origin as that of the short-ranged case: the indication is that second order perturbation 
theory is qualitatively good even for r^ ~ 1. 
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FIG. 17. (Ae)/A against Uo/t after 300 to 1000 disorder realisations, for a range of sample sizes, 
all at approximately quarter filling with W = 4. r^ ~ Uo/2t. The dashed line is proportional to 



Uq. 



To summarise the results for the mean Coulomb gap, we find that Koopmans' approx- 



imation ([Tl| ) makes an error in the mean charging energy which for small r^ and L, and 
fixed disorder scales like (Ae) oc r^L. There is also evidence that for the larger sizes 
{A > 50), far beyond that accessible by exact diagonalisation, that the size dependence 
vanishes: (Ae) oc r^. In contrast to the naive expectation however, we find no sign of this 
error vanishing relative to A in the thermodynamic limit. This is due to the fundamental dif- 
ference between occupied and unoccupied SCHF levels already discussed in the short-ranged 



case. 



We now consider the fiuctuations in A2. As an example of the interaction dependence 
of these fiuctuations in the various approximation schemes, we plot the results for a fixed 
size in figure [TB|. It is seen that applying Koopmans' theorem in the forms (|TI]-|T3p results in 
considerably smaller fiuctuations than the fully self-consistent calculation. To quantify this 
error, we plot (5A2^/(5A2 in Fig.|l9|, which shows that the relative error initially increases with 
interaction strength, but shows signs of saturating. The value of the saturation appears to 
increase towards unity as the system size is increased. 
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FIG. 18. Typical results for 5A2/A under various approximation schemes, averaged over 300 

disorder realisations. Here, W = 4, the lattice is 11 * 10 and N = 28, r^ ~ 0.54C/o/t. The 

Koopmans' approximants can be seen to underestimate the fluctuations. 
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FIG. 19. 6A2'^ /6A2 against Uq for a range of sample sizes, at approximately quarter filling: 
rg ~ Uo/2t. The statistics are obtained from 300 to 1000 disorder realisations for each sample size. 



We now concentrate on the fluctuations of the fully self-consistent peak spacing A2. For 
comparison with Ref. |^ it is useful to plot 8/^2/ {^2) against the interaction strength for a 
range of sample sizes. This is done in figure ^. In the inset we plot (5A2/A which shows 
that the peak spacing fluctuations are not proportional to A for VgL ~ 0{1) [^]. From Fig. 
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20|it can be seen that the curves 8/S.2/ {^2) do not saturate to a constant as suggested in Ref. 
0, although to see this clearly one has to consider larger sample sizes than are accessible 
by exact calculations. The curve 8/^2/ {^2) appears to take on the approximate form of a 
constant term plus a linear term for r^L ~ 0{1). The constant contribution identified by 
Sivan et el. 0], is here, contrary to their claim, non- universal (i.e. it is disorder dependent). 
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FIG. 20. (5A2/(A2) against U{)/t, after 300 to 1000 disorder realisations. The legend shows 
the sample size and W , v ~ 1/4: r^ ~ f7o/2t. The random matrix theory result at [/q = is 
8IS.2/ {IS.2) ~ 0.52. Inset: (5A2/A for the same data set. 



In figure ^ we plot 5^2/{^2) against disorder for the 10 x 11 lattice with a range 
of interaction strengths. At t/o = it is seen that for ly < 6 the systems obeys WD 
statistics quite well. For the sample size considered we find that in the regime 0.5 ^ 
f/o ~ 6.0 (0.25 ^ Ts ^ 3.0) 5A2/(A2) oc W , and at stronger interactions this dependence 
weakens. The intermediate dependence, 5A2/ {/S.2) oc W ., is consistent with the dependence 



6^2/ {^2) oc 1/y^ recently observed independently by Bonci and Berkovits ||30| for the 
Buminovich Stadium billiard. Analysis of Fig. ^ leads to the conclusion that the quadratic 
contribution (in Uq) to 6A2 is independent of disorder, which is consistent with Fig. |2^. 
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FIG. 21. JA2/(A2) against W/t averaged over 300 10 * 11 samples with N = 28. rs ^ O.hAUo/t. 

Results for other sample sizes were similar. The CI + RMT result for C/q = is plotted as a dotted 

line. 
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FIG. 22. 5A2/A against sample area A averaged over 300 to 1000 disorder configurations. The 
dotted line shows the CI + RMT result, the dashed line is proportional to yA, and the dot-dashed 
line is proportional to A. r^ ~ Uo/2t. 

To identify the system size scaling of the various contributions we plot 6A2/A against 
A in figure ^. For Uq = the system obeys WD statistics and 6A2/A is independent of 
size. The regime over which the fluctuations are approximately proportional to the mean 
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charging energy (the constant contribution to 8/^2/ {^2) alluded to above, which corresponds 
to a yA dependence in the figure), depends on the system size. As could be seen in Fig. 



20| another term begins to dominate the the fiuctuations at larger Uq which is quadratic in 
t/o, this term increases more rapidly with the system size, and so dominates at lower Uq in 
larger systems. Over the range of sizes considered, this term appears to scale like L^ . The 
cross-over in dominance therefore occurs at Uq ~ 1/L for fixed disorder strength; clearly 
the quadratic term will dominate in large samples. The increase in (5A2/A with system size 
appears to be an artifact of using the unscreened Coulomb interaction in the Hamiltonian. 



Summarising the results presented in figures ^ to ^, and the above discussion, we 
find for TgL > 0{l) an approximate form: 5/S.2 ~ .52A + a{/S.2) / \/g + br^ where a, 6 
are constants. One would normally expect the fiuctuations to be linear in the interaction 
strength, (i.e. 6 = 0). A possible source for such a quadratic interaction dependence in 
the typical fiuctuations is the development of correlations that grow like r^ in products of 
eight wavef unctions. Elsewhere [^ we present evidence for increased fiuctuations in the 
ground state density in this regime, as compared to a non-interacting system. It is not 
yet clear whether this result is an artifact of the SCHF approximation (which has also 
very recently been observed in Id systems PT| using a similar approximation scheme), or a 



genuine physical effect. 

Finally, it is worth noting that since the exchange interaction is not correctly screened, 
that errors in the SCHF scheme might be expected to diverge with respect to A as the system 
size is increased. We can neither confirm nor counter this argument, but have verified that 
for r^ ^ 5 the fiuctuations in the exchange contribution are smaller than those of the direct 
contribution. 

V. SUMMARY 

We have investigated the addition spectra of disordered quantum dots employing an 
effective single particle approximation, both using a fully self consistent analysis, and by 



30 



invoking Koopmans' theorem. We were able to consider system sizes with up to 144 sites, 
and 37 particles, compared to the latest exact calculations on samples with 24 sites and 6 
particles . The larger sample size also allows us to consider smaller values of g than exact 
calculations whilst retaining an ergodic non-interacting limit, and therefore approaches the 
experimental parameters more closely. The inclusion of spin in a consistent manner is left 
for a future project. 

Our SCHF results for the typical fluctuations of the peak spacings for particles possessing 
short-ranged bare interactions are entirely different from the results for long-ranged bare 



interactions. In the short-ranged case we find the same scaling as Ref. |I6| for very weak 
interactions (r^ -C 1), but for r^ ^ 1 deviations from this behaviour become significant and 
coincide with the appearance of interaction induced density modulations. We find no size 
dependence in the onset of these effects. We find that strong filling factor and geometry 
dependences arise due to these density fluctuations, and therefore do not expect that the 
disorder ensemble statistics can be mapped to statistics over the ensemble of filling factors: 
ergodicity is lost. We suggest that employing a short ranged bare interaction in a self- 
consistent scheme is not an appropriate model for the quantum dots of Refs. |^-^,|TT]-|I1| for 
which r^ > 1, but may be a useful model for dot geometries sandwiched between very close 
metallic gates which provide a good external source of screening [^,^. In this respect we 
identify some experiments on the addition spectrum which possess a metallic source (heavily 
doped n^ GaAs) and drain (Cr/Au), at separations of the order of the average inter-particle 
spacing in the dot [|,0. 

In the Coulomb case the SCHF approximation to A2 yields typical fluctuations that do 
not scale with A for r^ ^ 1/L. In Ref. it is claimed that 5/S.2 is universally proportional 
to (A2) for strong interactions (but still far from the accepted Wigner Crystal transition 
point). In contrast, we find, in addition to the small interaction independent contribution, 
a contribution to 8/^2 that is proportional to (A2)/y^ (i.e. non- universal) , and a further 
contribution that scales like r^, which is independent of disorder, and appears to be due 
to the development of charge density modulations |2^. The latter is not detectable in the 
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small systems examined numerically in |Q, and so our results are not numerically inconsis- 
tent with exact calculations [Q]. Whilst we do not include spin, the observed decrease in 
the fluctuations with g is consistent with the experimental indications [Q,^ that in cleaner 
samples the fluctuations are smaller. 

We show that a direct application of Koopmans' theorem overestimates (A2). This 
overestimate, a manifestation of the breakdown of Koopmans' approximation, does not 
vanish on the scale of A in the thermodynamic limit. The error seems to scale differently 
with sample size for sample areas above or below A k, 50. In the nearest neighbour case, 
with A ^ 50, this error scales with A, but in smaller systems, accessible by exact methods, 
it is independent of system size. In the Coulomb case the error grows with the system 
size as L for A ^ 50. For larger sizes the error appears to tend towards a 1/L scaling, 
i.e. in proportion with the charging energy, and therefore diverges with respect to the 
mean effective single-particle level spacing. This result for the Coulomb interaction case 
appears to be non-physical, and may be an artifact of the model considered. However, the 
result that Koopmans' theorem is not recovered in the thermodynamic limit also occurs 
in the short-ranged interaction case. In both cases we flnd that initially this error grows 
in proportion to f/g, to be expected since the lowest order contribution is second order, 
but for strong interactions it grows more slowly in [/q, and that the disorder dependence 
of this error is weak and non-monotonic. We identify the source of the error (Ae) to be 
the fundamental difference between occupied and unoccupied states that is inherent in the 
SCHF approximation. We introduce two improved applications of Koopmans' theorem, 
(Ag^), (Ag^), which provide a good approximation to (A2), but not to (5A2). 

Whilst preparing the manuscript, two related works appeared that conflrm some of the 
points discussed above [pn,Rl]. 



In both cases, fluctuations in the ground state density develop with r^ |Q, and have 
signiflcant effects of the addition spectrum statistics. It remains to be seen whether these 
density modulations are an artifact of the SCHF approximation (i.e. due the neglect of dy- 
namical correlations), or in fact interesting results on the continuous transition to a Wigner- 
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type solid in disordered samples with short- and long-ranged bare interactions. 
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